home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / kaos / kaos_pat.0 / trkman_ode.c < prev   
C/C++ Source or Header  |  1991-10-04  |  3KB  |  118 lines
  1. /*
  2. ----------------------------------------------------------------------
  3.      Compute a set of manifolds for saddle equilibria of ODE
  4. ----------------------------------------------------------------------
  5.      compute manifolds of ANY given set of hyperbolic fixed points
  6.      with eigenvalues and eigenvectors
  7.      iskip : number of skips between each data
  8. Input:    sevec[n][n]: eigenvectors for n eigenvalues
  9.     seval[n]: eigenvalues of an equilibrium
  10.     xe[n]: coords of an equilibrium
  11.     mu: # of maximum steps along unstable manifold
  12.     ms: # of maximum steps along stable manifold
  13.     muf : # of finer division of the unstable manifold
  14.     msf : # of finer division of the stable manifold
  15.     iskip: # of skips before drawing
  16.     delman: initial distance from an equilibrium 
  17.     n: phase space dimension
  18.     r: unused, always set to 0 within the subroutine (kept to retain the
  19.     same syntax as trkman_map.c)
  20. ----------------------------------------------------------------------
  21. NOTE:    -Only works for r=0; (equilibrium point of ODE)
  22.     -Internal variable: man_index: -1=complex, 0=unstable, 1=stable
  23. ----------------------------------------------------------------------
  24. */
  25. #include "../include/x11r2_kaos_def.h"
  26.  
  27. #define RK4 1    /* Runge-Kutta integrator */
  28. #define DELFRAC 0.1
  29. void trkman_ode(sevec,seval,xe,r,ms,mu,msf,muf,iskip,delman,n)
  30. int r,ms,mu,msf,muf,iskip,n;
  31. double **sevec,**seval,*xe,delman;
  32. {
  33.     int i,ic,jc,mc,mend,man_index,icnt,color;
  34.     double fabs(),t_step,delx,*x,*xp,*x1,*x2,*dvector();
  35.     extern int stop,int_algorithm,dot_size;
  36.     extern double time,time_step;
  37.     extern char string[];
  38.  
  39.     r=0; /* r should be 0 */
  40.     x1 = dvector(0,n-1);
  41.     x2 = dvector(0,n-1);
  42.     x = dvector(0,n-1);
  43.     xp = dvector(0,n-1);
  44.  
  45.     for(ic=0;ic<n;ic++){    
  46.         if(seval[ic][1]!=0) /* if complex eigenvalue quit*/
  47.             man_index = -1;
  48.         else {
  49.             if(seval[ic][0]>0){
  50.                 man_index = 0;
  51.                 color = Blue;
  52.             }
  53.             else {
  54.                 man_index = 1;
  55.                 color = Red;
  56.             }
  57.         }
  58.         if(man_index== 0){
  59.             t_step = time_step/muf; /* temporary time step */
  60.             mend = mu * muf;
  61.         }
  62.         else if(man_index == 1){
  63.             t_step = - time_step/msf;
  64.             mend = ms * msf;
  65.         }
  66.  
  67.         sprintf(string,"Computing manifold for ev[%d]...",ic);
  68.         printf("%s\n",string);
  69.         printf("(Ev=%g %g)\n",seval[ic][0],seval[ic][1]);
  70.         printf(" Evec %d = %g \n",ic,sevec[ic][0]);
  71.         for(i=1; i<n; i++)
  72.         printf("          %g \n",sevec[ic][i]);
  73.  
  74.         for(jc=0;jc<2;jc++){
  75.             if(man_index==0){
  76.                 printf("Unstable manfold for ev[%d]...\n",ic);
  77.             }
  78.             else if(man_index==1){
  79.                 printf("Stable manfold for ev[%d]...\n",ic);
  80.             }
  81.             if(jc==0){
  82.                 for(i=0;i<n;i++) x[i] = xe[i]+delman*sevec[ic][i];
  83.             }
  84.             else {
  85.                 for(i=0;i<n;i++) x[i] = xe[i]-delman*sevec[ic][i];
  86.             }
  87.             icnt=0;
  88.             for(mc=0;mc<mend;mc++){
  89.                 for(i=0;i<n;i++){
  90.                     delx=(x[i]-xp[i])/seval[ic][0];
  91.                     x2[i]=x[i]+delx;
  92.                     x1[i]=x1[i]- DELFRAC*delx;
  93.                 }
  94.                 /*use Runge-Kutta 4th order */
  95.  
  96.                 int_onestep(x1,x2,x,&time,t_step,n,RK4);
  97.                 for(i=0;i<n;i++){
  98.                     xp[i] = x[i];
  99.                     x[i] = x1[i];
  100.                 }
  101.  
  102.                 if((icnt % iskip)==0){
  103.                     (void) draw_record_pwf(x1,color,-1,dot_size,0,1);
  104.                 }
  105.                 icnt++;
  106.                 if(stop){
  107.                     goto done;
  108.                 }
  109.             }
  110.         }
  111.     }
  112. done:
  113.     (void) free_dvector(x1,0,n-1);
  114.     (void) free_dvector(x2,0,n-1);
  115.     (void) free_dvector(x,0,n-1);
  116.     (void) free_dvector(xp,0,n-1);
  117. }
  118.